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^ Abstract 

CSJ We study the dynamics of a pair of parametrically-driven coupled nonlinear mechanical res- 

(N 

onators of the kind that is typically encountered in applications involving microelectromechanical 

\-A and nanoelectromechanical systems (MEMS & NEMS). We take advantage of the weak damping 

_H that characterizes these systems to perform a multiple-scales analysis and obtain amplitude equa- 

'^ tions, describing the slow dynamics of the system. This picture allows us to expose the existence 

of homoclinic orbits in the dynamics of the integrable part of the slow equations of motion. Using 

^ a version of the high-dimensional Melnikov approach, developed by Kovacic and Wiggins [Physica 

^I^ D, 57, 185 (1992)], we are able to obtain explicit parameter values for which these orbits persist 

in the full system, consisting of both Hamiltonian and non-Hamiltonian perturbations, to form so- 
^^ called Silnikov orbits, indicating a loss of integrability and the existence of chaos. Our analytical 

^\^, calculations of Silnikov orbits are confirmed numerically. 
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I. INTRODUCTION 

Microelectromechanical system (MEMS) and nanoelectromechanical systems (NEMS) 
have been attracting much attention in recent years [IHS]. MEMS & NEMS resonators are 
typically characterized by very high frequencies, extremely small masses, and weak damping. 
As such, they are naturally being developed for a variety of applications such as sensing with 
unprecedented accuracy |1H6], and also for studying fundamental physics at small scales — 
exploring mesoscopic phenomena [71 |8] and even approaching quantum behavior [9l [TO] . 
MEMS & NEMS resonators often exhibit nonlinear behavior in their dynamics [II1II2]- This 
includes nonlinear resonant response showing frequency pulling, multistability, and hystere- 
sis [31 fT3l - [T6] . as well as the formation of extended [17] and locahzed [18] collective states 
in arrays of coupled nonlinear resonators, and the appearance of chaotic dynamics [T9H2T] . 
Nonlinearities are often a nuisance in actual applications, and schemes are being developed 
to avoid them [22], but one can also benefit from the existence of nonlinearity, for example 
in mass-sensing applications [231 [21], in achieving self-synchronization of large arrays [23], 
and even in the observation of quantum behavior [2B] • 

MEMS & NEMS offer a wonderful experimental testing ground for theories of chaotic 
dynamics. Numerical investigations of a number of models of MEMS & NEMS resonators 
have demonstrated period-doubling transitions to chaos [211 [27[430] , yet there are very few 
analytical results. One of the simplest models of chaotic motion is that of the Duffing 
resonator with a double-well potential, described by the Hamiltonian 

with k < and a > 0. This simple mechanical system has a homoclinic orbit for H = 0, 
connecting a saddle at the origin of phase space to itself. Upon the addition of damping and 
an external drive this system develops a particular kind of chaotic motion called horseshoe 
chaos [31], which can be studied analytically using the Melnikov approach [32] • The stable 
manifold leading into the saddle and the unstable manifold leading away from the saddle, 
which coincide in the unperturbed Hamiltonian ([I|, are deformed when damping and a drive 
are added. Yet, conditions can be found analytically using the Melnikov function, which 
measures the distance between the two manifolds, under which they intersect transversely 
leading to the possibility of observing chaotic dynamics. What one observes in practice is 
a random-like switching of the resonator between the two wells. Thus, having an analytical 



criterion for asserting the existence of chaotic motion allows one to distinguish it from 
random stochastic motion that might arise from noise. Such a Hamiltonian as in Eq. ([I]) 
was implemented in a MEMS device using an external electrostatic potential by DeMartini 
et al. [20] , and studied using the Melnikov approach. 

Here we wish to study the possibility of observing horseshoe chaos in typical NEMS 
resonators, which are described by a potential as in Eq. ([I]), but of an elastic origin with k 
and a both positive. Individual resonators of this type do not exhibit homoclinic orbits for 
any value of H, and therefore are not expected to display horseshoe chaos under a simple 
periodic drive. Nevertheless, a pair of coupled resonators of this kind — like the ones studied 
experimentally by Karabalin et al. [21] — are shown below to possess homoclinic orbits in 
their collective dynamics, and are therefore amenable to analysis based on a high-dimensional 
version of the Melnikov approach [33|. We employ here a particular method, developed by 
Kovacic and Wiggins [31| , which is a combination of the high-dimensional Melnikov approach 
and geometric singular perturbation theory. This method enables us to find conditions, in 
terms of the actual physical parameters of the resonators, for the existence of an orbit in 4- 
dimensional phase space, which is homoclinic to a fixed point of a saddle-focus type. Such an 
orbit, called a Silnikov orbit [35j, provides a mechanism for producing chaotic dynamics [31]. 

We study here the case of parametric, rather than direct driving, but this is not an 
essential requirement of our analysis. On the other hand, having weak damping, or a large 
quality factor, characteristic of typical MEMS & NEMS resonators, is essential for the 
analysis that follows. First of all, as was demonstrated in a number of earlier examples [36l 
ET] , it leads to a clear separation of time scales — a fast scale defined by the high oscillation 
frequencies of the resonators, and a slow scale defined by the damping rate. This allows us 
to perform a multiple-scales analysis in Sec. |TT] and obtain amplitude equations to describe 
the slow dynamics of the system of coupled resonators. It is in the slow dynamics that the 
homoclinic orbits are found. Secondly, the weak damping, which requires only a weak drive 
to obtain a response, allows us to treat both the damping and the drive as perturbations. 



even with respect to the slow dynamics. Therefore, in Sec. Ill we set the parametric drive 



amplitude and the damping to zero in the amplitude equations, which makes them integrable. 



This allows us, in Sec. |IV[ to find conditions for the existence of homoclinic orbits and 
to obtain analytical expressions for these orbits. We emphasize that these orbits reside 
in a 4-dimensional phase space, and as such are homoclinic not to a point, but rather 



to a whole invariant 2- dimensional manifold in the shape of a semi-infinite cylinder. Of 
these, we identify a subset of orbits, satisfying a particular resonance condition, that are 
precisely heteroclinic, connecting pairs of points in 4-dimensional phase space. In Sec. |V] 
we reintroduce the drive and the damping into the equations as small perturbations, and 
use the high- dimensional Melnikov method to determine which of the heteroclinic orbits, 
determined through the resonance condition in the unperturbed system, survives under the 



perturbation. In Sec. |VI| we study the effects of the perturbation on the dynamics within 

we 



the invariant semi-infinite cylinder near the resonance condition. Finally, in Sec. VII 
put everything together by calculating the parameter values for which the end-points of 
the unperturbed heteroclinic orbits are deformed in the perturbed system in such a way 
that they become connected through the dynamics on the semi-infinite cylinder, producing 
Silnikov orbits, homoclinic to a fixed point of a saddle-focus type. We conclude by verifying 
our analytical calculation using numerical simulations. Our analysis implies that conditions 
exist in the coupled resonator system that could lead to chaotic motion. 

II. NORMAL MODE AMPLITUDE EQUATIONS 

We consider a pair of resonators modeled by the equations of motion 

Un + Un + ul- ^Q'^{Un-l " 2m„ + Un+l) 

+ -[D + HcOSUJpt]{Un-l -2Un + Un+l) = 0, n = 1, 2, (2) 

where m„ describes the deviation of the n^^ resonator from its equilibrium, and we label 
two fictitious fixed resonators as uq = M3 = for convenience. Detailed arguments for 
the choice of terms introduced into these equations of motion are discussed by Lifshitz and 
Cross [36], who modeled the particular experimental realization of Buks and Roukes [T7] . 
although other variations are possible [TT]. The terms include an elastic restoring force as in 
Eq. ([I| with positive linear and cubic contributions (whose coefficients are both scaled to 1), 
a dc electrostatic nearest-neighbor coupling term with a small ac component responsible for 
the parametric excitation (with coefficients D and H respectively), and a linear dissipation 
term, which is taken to be of a nearest neighbor form, motivated by the experimental 
indication [TTj that most of the dissipation comes from the electrostatic interaction between 
neighboring beams. Note that the electrostatic attractive force acting between neighboring 



beams decays with the distance between them, and thus acts to shghtly soften the otherwise 
positive elastic restoring force. Lifshitz and Cross also considered an additional nonlinear 
damping term, which we neglect here for the sake of simplicity. The resonators' quality 
factor Q is typically high in MEMS & NEMS devices, which can be used to define a small 
expansion parameter e <^ 1, by taking Q^^ = e^, with 7 of order unity. The drive amplitude 
is then expressed as if = eh, in anticipation of the fact that parametric oscillations at half 
the driving frequency require a driving amplitude which is of the same order as the linear 
damping rate [TT]. 

Following Lifshitz and Cross |36l Appendix B] , we use multiple time scales to express the 
displacements of the resonators as 

xi,2{t) = ^ (Ai(T)e^-^* ± A2{T)e'^'' + c.c.) + e^/'xg(t) + ..., (3) 

where Xi is taken with the positive sign and X2 with the negative sign; with a slow time T = 
et, and where the normal mode frequencies are given by w^ = 1 — D/2, and ijj\ = \ — 'iD jl. 
Substituting Eq. (|3]) into the equations of motion (|2| generates secular terms that yield 
two coupled equations for the complex amplitudes A\^2- If we measure the drive frequency 
relative to twice 1^2 by setting uj^ = 2uj2 + efi, express Ui relative to UJ2 as uji = uj2 + 2ef2i, 
and express the complex amplitudes using real amplitudes and phases as 

the real and imaginary parts of the two secular amplitude equations become 

ddi 1 ^ h 9 2 / \ /■<- \ 

-^ = -77^1 - o — «ism2xi - - — a2«ism2(x2 - Xi), (5a) 

CLl 4 OUJi oWi 

^ = 2^1 - ^fi - -^ cos 2x1 + ^ K + 2^2 + «2 cos 2(x2 - Xi)] , (5b) 

da2 3^ 3h ^ 2 ■ of ^ ^k ^ 

-p^ = -7702- ^a2sm2x2- ^0^02 sm2(xi - X2), (5c) 

CLl 4 0UJ2 0W2 

"X.2 ^ f~\ ^'^ 9r2 22/ \1 /i\ 

-7^ = --^- ^cos2x2 + ^[a2 + 2ai + aiCos2(xi-X2)J- (5d) 

CLl Z 0UJ2 0UJ2 

Steady-state solutions, oscillating at half the parametric drive frequency, are obtained by 
setting dai/dT = dxi/dT = in Eqs. (|5| and solving the resulting algebraic equations. We 
are interested in extending the investigation of these amplitude equations. In particular, 



we want to identify the conditions under which they may display chaotic dynamics. We 
should note that equations similar to (tsl) were also used for modeling a variety of paramet- 
rically driven two-degree of freedom systems such as surface waves in nearly-square tanks or 
vibrations of nearly-square thin plates or of beams with nearly-square cross sections 



III. UNPERTURBED EQUATIONS SETTING DAMPING AND DRIVE TO 
ZERO 



We first consider the integrable parts of Eqs. (pj), obtained by setting 7 = /i = 0, which 
after a rescaling of the amplitudes, ai —^ ai ^yu28/9, and 02 — )■ 02 \/ui8/9, become 

— ^ = alai sin 2(xi - X2), (6a) 

'^^'- ^^ - 2n,) + a^ [2 + cos 2(xi - X2)] + -a?, (6b) 



ala2sm2{xi-X2), (6c) 



(i(i2 2 



dT 

^ = -!^ + «? [2 + cos2(xi - X2)] + ^al (6d) 

aJ z (jj2 

In Sec. |V] we will reintroduce the driving and damping terms as a perturbation. We trans- 
form Eqs. (Ig]) into a more familiar form, which has been studied in the context of higher 
dimensional Melnikov methods ^ST[ El], by changing to two pairs of action-angle variables: 
(i) B = al/2, = xi — X2', and (ii) I = {af + a\)/2, (p = X2- After defining 6 = uJi/u2, and 
rescaling time as T — ;■ T/2, we obtain the unperturbed Hamilton equations 

§ = ^^'f^^'^^ =n, + I{2-6 + cos2e)-B(^A-^-^ + 2cos2ey (7b) 

dl _ dHo{B,e,I) 
dT~ dd) 



(7c) 



d0 dHoiB,9,I) O 

— = — = SI-- + B{2-S + cos2e), (7d) 

where the Hamiltonian Hq, which generates these equations, is expressed as 

H,{B,9,I) = ^-^-^I- B' (2 - ^^"j + B [/(2 -6) + n,]+B{I- B) cos 20. (8) 

Thus, both I and Hq are constants of the motion in the unperturbed system. Note that 
[B, 6, 1, 0) G M^ X S X ]R+ X §, where S is the unit circle, and IR+ are the non-negative reals. 



It is convenient to describe the dynamics also in terms of the Cartesian variables x = 
ai cos(xi — X2) = V'2B cos 9 and y = ai sin(xi — X2) = V2Bsm9, in place of B and 9, 
thereby obtaining the Hamilton equations 

dx dHo{x,y,I) 3/ 6^ + l\ ^ fn ^^ + 1 



y^(l-^)+^^y(2-^)-y[/(l-5) + f^i], (9a) 
'' - '^°^^' ^' '^ - x^ fa - ^ V y^x ^2 - ^ V ^ [/(3 - S) + ..] , (9b) 



dT dy 



dT dx V 2(5 / ^ V 2(5 

dl dHo{x,y,I) 



dT d(f) 



0, (9c) 



^-^^^^-^/-f + y(3-^)4(l-^), (9d) 

where y plays the role of a coordinate and x is its conjugate momentum, and where the 
Hamiltonian Hq is now given by 

+ y m - '5) + ^1] + Y m -s)+ ^1] • (10) 



IV. ANALYTICAL EXPRESSIONS FOR HOMOCLINIC ORBITS 

We wish to identify the conditions under which there exist homoclinic orbits in the unper- 
turbed system. These orbits will potentially lead to chaotic dynamics once we reintroduce 
the damping and the drive in the form of small perturbations. We therefore consider the 



fixed point x = y = in the unperturbed {x,y) plane, as given by Eqs. (9a) and (9b). A 
linear analysis of this fixed point reveals that it is a saddle for values of the positive constant 
of motion J, that satisfy the inequality [/(I — (5) + fii][J(3 — 6) + Qi] < 0. This implies that 
the fixed point at x = y = is never a saddle if the fixed parameter 5 < 1; it is a saddle for 
1 < (5 < 3, if / > fii/((5 - 1); and it is a saddle for (5 > 3, if fii/(5 -!)</< fii/(5 - 3). 
We shall restrict ourselves here to values 1 < (5 < 3, therefore to obtain a saddle one must 
only ensure that / > f2i/(5 — 1). In the full four-dimensional system given by Eqs. ^ this 
saddle point describes a two-dimensional invariant semi-infinite cylinder, or annulus, 

(x,y,J,0) \x = 0, y = 0, ^</}, 1<(5<3, (11) 

where is unrestricted within the unit circle. The trajectories on ^ are periodic orbits 
given by / = constant and = {61 — Q/A)T + 0o- For the resonant value of / = f = Q/AS 



the rotation frequency vanishes, and the periodic orbit becomes a circle of fixed points. Of 
course, this trivial unperturbed dynamics on ^ undergoes a dramatic change under the 
addition of perturbations. 

The two-dimensional invariant annulus ^ has three-dimensional stable and unstable 
manifolds, denoted as W''^{^) and W^{^), respectively, which coincide to form a three- 
dimensional homoclinic manifold F = W'^{^) fl W^{^). Trajectories on the homoclinic 
manifold F are homoclinic orbits that connect the origin of the (x, y) plane to itself. Thus, 
the constant value of the Hamiltonian along such an orbit is equal to its value at the origin, 
namely Hq{Q, 0, /) = Hq{x, y, I) = Hq{B, 9, 1), which immediately yields an equation for the 
homoclinic orbits in terms of the action-angle variables 

^ ^^' ^^ - 52-25(2 + cos20) + l ■ ^^^' 

To obtain the temporal dependence of the dynamical variables along the homoclinic orbit. 



we substitute the homoclinic orbit equation (12) into Eq. (7b), to get 

de 



I{6 -2 -cos2e)-ni. (13) 



Next, we note that Xi = + ^; and use the Hamiltonian dsj) to get 



f--?-^- (") 



We then integrate Eq. (13), substitute the result into Eq. (12), and the latter into Eq. (14) 



and finally integrate Eq. ( 14 ) to obtain analytical expressions for the temporal dependence 
of the dynamical variables along orbits that are homoclinic to ^. 

For / > 25ni/{6^ - 1) we define q = I{5^ - 1) - 25^1 > 0, and find that Oq = 9{T = 0) = 
0, IT, and that the homoclinic orbits are given by 

B\T, I) = -^^^— , (15a) 

^ ' ' gcosh(2aT)+p' ^ ' 

tan {e'^iTJ)) = -y f(iI^|~[^| tanh(aT), (15b) 

X1{T,I) = -44= arctanh (" /^tanhar) + (6 1 - ^)T + xi{0) , (15c) 
<P\T,I) = x1iT,I)-e\T,I), (15d) 

p = Q^{S^ _ 45 + 1) _ 7(^3 _ Q^2 ^ 7^ _ 2)^ 
a^ = -nl + 2/fii((5 - 2) - I\6 - 3){6 - 1). (16) 



where 




(a) / = 6 
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FIG. 1. Numerical phase portraits of the unperturbed (x,y) plane for 5 = 2, fii = 21.32, and 
different values of I as noted in the individual captions. All figures show the trajectories for which 
B has a fixed value equal to I. As explained in the text and for 5 = 2, four hyperbolic fixed points 
appear on the B = I circle for 2J7i/5 = 8.528 < I < 2ili = 42.64; the B = fixed point changes 
from a center to a saddle at / = ili = 21.32; and a global bifurcation rotating the homoclinic orbit 
through an angle of it/2, shown in panel (d), occurs at / = 4r2i/3 ~ 28.4267. 

For / < 26ni/{6^ - 1) we redefine q = 26Qi - I{6^ - 1) > 0, and find that 9o = ±n/2, 



(17) 



and that the homoclinic orbits are given by Eqs. (Il5j), with Eq. (15b) replaced by 

Thus, exactly at / = 26Qi/{6'^ — 1) (or q = 0) there is a global bifurcation in which the 
homoclinic orbit rotates through an angle of tt/2. 

Some of the phase-space portraits of the unperturbed (x, y) plane are calculated numer- 
ically from Eqs. (7a) and (7b), for different values of /, and shown in Fig. [I] From Eq. (7a) 
it follows that the value of B is fixed if (a) B = 0; or (b) B = I; ot (c) 6 is an inte- 
ger multiple of 7r/2 and 6 is fixed. Figures [ita)-(f) all show the trajectories for which B 




FIG. 2. Orbits homoclinic to ^ . For / = V (the orbit in the middle), d(j)/dT = on ^, and the 
orbit is heterochnic, connecting fixed points on ^ that are A.(j) apart. For / ^ r\ d4>/dT ^ on 
^ . The parameters are 5 = 2,Vt = 400, i7i = 21.32. 



has a fixed value equal to /. Four hyperbolic fixed points appear on the B = I circle for 
r2i/(3 — 1/5) < I < fii(l — 1/5), where solutions exist to the equation dO/dT = with 
B replaced by / [Figs. [I](b)-(e)]. As expected, the origin i? = is always a fixed point — a 
center for small values of I [Figs. [l](a),(b)], which undergoes a pitchfork bifurcation into a 
saddle when 89 /dT = with B = 0, occurring at / = rii/{5 — 1) [Figs. [T](c)-(f)]. Additional 
centers appear whenever 6 is an integer multiple of 7r/2 and solutions exist to the equation 
89 /dT = with cos 2^ = ±1 [Figs. |lfb)-(f)]. The global bifurcation at / = 2^15 /{6^ - 1) 
where the homoclinic orbit rotates by 7r/2 is shown in Fig. flld). 

Note that we refer to the orbits given by Eqs. (fTsl) as homoclinic since they are homoclinic 
to ^. A few of these orbits are shown in Fig. |2} At resonance, for / = f, the orbits are 
truly heterochnic, connecting fixed points that are A0 apart, where A0 = Axi — A6, and 



A^ = — 2arctan4 



Pil-6) + Qi 



Axi 



2a(52-l) ip-q 

-arctanhi 



(18a) 
:i8b) 



and where for any variable /, A/ = f{T = oo) — /(T = — oo). Such an unperturbed 



10 



10 




a^sinlx^-Xg! 



'=;::::b 




FIG. 3. Results of a numerical integration of Eqs. (pi) for different values of the initial amplitude 
of the second mode 02(0), with A = 7 = 0, ^i = 0.8528, uj2 = uJi/2, Q.i = 21.32, and Q. = 400. As 
expected, for 02(0) > Y^16riia;ia;2/9(a;i —(^2) — 5.68 the origin becomes a saddle, which rotates 
through an angle of 7r/2 at 02(0) = -y/32$7ia;fa;2/9(a;J — a;|) ~ 6.46. 



heteroclinic orbit is shown in the middle of Fig. |2} 



We wish to demonstrate the results obtained so far also in terms of original amplitude 
equations ([5| with /i = 7 = 0. The point x = y = corresponds to ai = in Eqs. (|5J). 
To start the simulation near this point, we initiate the numerical solution with ai(0) ^ 
1, which through the definition of / implies that 02(0) ^ yl6/cJi/9. The condition for 



having a saddle at the origin of Eqs. ( |9a| and (9b), / > ^i/{5 — 1), translates into the 
condition 02(0) > ^/lQVl^uJ]J^{5~^^V) = A/l6f2ia;iW2/9(a;i — 102)- The condition for the 
global bifurcation, rotating the homoclinic orbit through 7r/2, given by / = 25VLi/{5'^ — 1), 



translates into 02(0) = A/32r2iu;^u;2/9(a;f — a;|). These conditions are verified by a numerical 
integration of Eqs. ([s]) by varying the initial amplitude of the out-of-phase mode, 02(0), as 
shown in Fig. [3j 



11 



V. HOMOCLINIC INTERSECTIONS IN THE PERTURBED SYSTEM 

After having calculated the homoclinic orbits in the unperturbed system, we now reintro- 
duce the drive and the damping as perturbations and study how they affect the dynamics. 
In particular, we want to study the nature of the invariant annulus ^, and its stable and 
unstable manifolds, W'^{^) and W^"(^), under the perturbation, and use the Melnikov 
criterion to find the conditions under which they can still intersect. The perturbed equations 
are written in terms of the action-angle variables as 



— = 2B{I - B) sin 29 - ^hB sin 2(0 + 9)- ^-/B, (19a) 

-— = ni + 1(2 - 6 + cos29) - B { 4 ^ + 2 cos 29 

dT \ 

- -y [cos2(0 + 0)-35cos20], (19b) 

^ = -^h [35(/ - B) sin 20 + 5 sin 2(0 + 9)] - ^^{31 - 2B), (19c) 

^ = -^ + 5/ + 5(2-5 + cos2^)-^cos(20), (19d) 



where we have quantified the perturbations by expressing the drive amplitude and the damp- 
ing as h = C,Suih and 7 = ,^47, respectively, where ^ <^ 1 is a small parameter. It is 
instructive to write the perturbed system in the general form 



dB _ dHo{B,9,I) , ^ g_ dJUBAJl dH,{B,9,IA) , ,b\ ... ^ 

dT- d~9 +^^ " d~9 ^^[ d~9 ^"^ J' ^^"^^ 

d9 dHo{B,9,I) , dHoiB,9,I) ^ dH,iB,9,I,<P) 

dT = ^^ + ^^ =^ +^ dB ' ^2°^) 

rf0 dHoiB,9,I) ^ dHoiB,9,I) ^ dH,{B,9J,cp) 

dT = — dl — + ^^ - — dl — + ^ m ' ^'°^) 



where the perturbations due to the parametric drive are generated from the Hamiltonian 

iJi(5,^,J,0) = --[5cos2(0 + ^) + 35(J-5)cos20], (21) 

12 



and the dissipative perturbations are given by d^ = —^B and d^ = —7(3/ — 2B). Similarly, 
in terms of the Cartesian variables, the perturbed system is written in this general form as 
dx dHo{x,y,I) ^ dHo{x,y,I) ^ / dHi{x,yJ,4>) ^ \ 



dT dy dy \ dy 

dy dHo{x,y,I) ^ dHo{x,y,I) ^ / dHi{x,y,I,^, ^ ^^ . 



is' = ( (^-'J^i^^^ + d') 



dT dx dx \ dx 

dl 

dT 



(22c) 



with 



d^ _ dHo{x,y,I) , ^ ^_ dHo{x,y,I) ^ ^ dHi{x,y,I,^) 

dT~ dl +^^ " dl +^ dl ' ^ ' 

Hi (x, y, /, 0) = - { [(35 - 1) x^ + (3(5 + 1) y2 - Q5I] cos 20 + 2xy sin 20} , (23) 

and where d^ = —'yx/2, d^ = —'~fy/2, and d^ = —7(3/ — x^ — y^). 

For < ^ ^ 1 the unperturbed invariant annulus ^, and its stable and unstable 
manifolds, W'^{^) and W^{^), persist as a locally invariant annulus ^^ with stable and 
unstable manifolds, W'i^/:) and Wi^^) ^311 EH ESI I12]. Due to the fact that we use 
parametric rather than direct excitation, the point x = y = remains a fixed point of the 



perturbed Eq. (22a) and (22b), so ^^ is defined just like ^ in Eq. (11). However, the 
term locally invariant means that trajectories with initial conditions on ^^ may leave it 
through its lower boundary at / = Qi/{6 — l). We want to find intersections of the manifolds 
W^{^^) and W^{^^), because such intersections may contain orbits that are homoclinic to 
^^. This is done by calculating the Melnikov integral, M{I,(f)Q), which is a measure of the 
distance between these manifolds. If the Melnikov integral has simple zeros [M{I, 0o) = 
and (9M(/, 0o)/90o 7^ 0], the three-dimensional manifolds W^{J^^) and W^{^^) intersect 
transversely along two-dimensional surfaces. 
The Melnikov integral is given by 



where 



/oo 
(n(x\ y\ I), g(x^ y\ J, 0^^ + 0o))rfT, (24) 

-00 



_ f dHoix,y,I) dHoix,y,I) dH,{x,y,I) dHoiO,0,I) \ 
'"^"''^'^^"V dx ' dy ' dl dl J' ^^^' 

g{x,y,I,<P) = {g^,gy,g'), (26) 



x'^{T,I), y^{T,I), and (j)^{T,I) are the homoclinic orbits given by Eqs. (15), and angular 



brackets denote the standard inner product. At resonance, the Melnikov integral M(/'',0o) 
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can be calculated explicitly, because then dHQ{0,0,I^)/dI = and the integrand of the 
Melnikov integral is given by 



dHo OHq dHo J 

^ _m,9Hi imdHi _ dihdHi <m^, dihj, ago J 

dx dy dy dx dl dcj) dx dy dl 



(27) 



For the unperturbed orbits we can use the chain rule and the fact that dl/dT = to 
obtain the relation 



dHi dHodHi dH^dH^ dH^dHi 



dT dx dy dy dx 
so the Melnikov integrand reduces to 



dl 



dH^ dH. dH. dH^i 

Upon transforming to the action-angle variables one has 



(28) 



(29) 



dx 



dy 



and so the integrand (29) becomes 

(n, g) = 



dH^ dHo^B^dHo^j 
+ -T^ — d H — TT^d 



dT dB 
dHi 



dl 



dB 



.^-7S--7(3r-25)^ 
dT ' dT '^ Ut 



(30) 



-37/^^ + 275^-375-, 
dT ' dT ' dT ' dT' 



(31) 



where we recall that Xi = ^ + 



We can now explicitly integrate each of the terms in the integrand (31). From Eq. (21), 



owing to the fact that on the homoclinic orbits i?(±oo) = 0, the first of these yields 



dHi 

It 



dT 



36rh 
36rh 



cos20(oo) — cos20(— 00)] 

cos 2(00 + A0/2) - cos 2(00 - A0/2)] 



361^ h sin 20o sin A0, 



(32) 



00 



-00). The second term in (31) immediately yields 



where we recall that A0 = 
-~3'~fI'^A(f). For the third term in (31) we use Eq. (14), which on resonance yields 

1-52 '"^ 



B^dT 
dT 



26 
[1 - 6^)26a^ 



B^dT 



2p 



(p — o^w/z 



arctanh^ / h 



p + q q^ — p^ 



Aa. (33) 
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For the fourth and last term in (31) we use Eq. (12) and get 

5-1 



I 



Bdd = re + 



r{6^ - 1) - 25^1 

{6 - l)VW^66Tl 



arctan 



V6^ -65 + 1 



tan^ 



(34) 



and after substituting the hmits, using Eq. (15b), we get 



rfx2 



Bd9 = FM - 2 



r{5 



25VLi 



= FAe + A/i. 



After collecting all four terms we finally obtain 



arctan 



'/'■(5-3)-fii 



VS^ -6(5+1 V ^'(1 - 5) + ^1 



(35) 



M{F,^o) = -35r/isin20osinA0-7(3rAxi + 3A^-2Aa). 



(36) 



Except for the special case in which the phase difference A0 is a multiple of vr, the function 
M{F, 0o) has simple zeros as long as the relation 



7(3r Axi + 3A/i - 2Aa) 



35Fh sin A0 



< 1 



(37) 



is satisfied. If the system parameters satisfy this condition, every simple zero of the Mel- 
nikov function corresponds to two symmetric (due to the invariance x,y ^ — x, —y) two- 
dimensional intersection surfaces. The ^ — )■ limit of these surfaces contain orbits whose 



explicit form is given by Eqs. (15), with their / and 0o values satisfying the relation 
M(J, 0o) = 0, for / close to /'' [M]. Thus, an unperturbed heteroclinic orbit given by 



Eqs. (15), with / = F and a phase 0o at time zero, can be made to persist under the 



perturbation by setting the drive amplitude to the value 

7(2Aa - 3FAxi - 3A/i) 



h 



36F sin 20o sin A0 



(38) 



We give numerical evidence of this in Sec. VII Such orbits surviving in the intersection 
of W^{^^) and W^{^^) may leave the stable manifold W^{^^) in forward time, and the 
unstable manifold W^{^^) in backward time, through the low boundary at / = fii/(5 — 1), 
since these manifolds are only locally invariant [34J. However, the analysis we perform 
below allows us to find surviving homoclinic orbits that are contained in the intersection of 
W{^^) and W 
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VI. DYNAMICS NEAR RESONANCE 

After having calculated the Melnikov integral at / = /'', we proceed to examine the 
dynamics on ^^ near this resonance. The equations that describe the dynamics on ^^ are 



obtained by setting i? = in Eqs. ( 19c ) and ( 19d ) 



dl 

df 

d(f) 

dT 



-^3/(Msin2(/) + 7) 
~4 



61 



M6 

4—— cos 20. 



To investigate the slow dynamics, which is induced by the perturbation on ^ 
nance, we follow Kovacic and Wiggins [34J and introduce a slow variable I = F 



Eq. (39), along with a slow time scale r = a/^T, and obtain 



^ = -3(r + V^p)(Msin20 + 7), 
or 

d(t> r /T3h5 



The leading terms in Eqs. (40), independent of ^, yield 



dp 
dcf) 



6p 



3r(Msin(20)+7) 
9^(p,0) 



d,jr{p, 



where 



^(p. 



dp 



]5 p^ -IhSr cos{2^) + S-fF 



(39a) 
(39b) 

near reso- 
f y/^p into 

(40a) 
(40b) 

(41a) 
(41b) 

(42) 



2 ^ 2 
is a rescaled Hamiltonian that governs the slow dynamics on ^^ close to resonance. 

Fig. |4]^a) shows the phase portrait of Eqs. (41 ), which contains a saddle go at (p = 0, = 
0s = [arcsinfe — 7r]/2), and a center Pq at (p = 0,0 = 0c = — [arcsin6]/2), where b = '~f/h6. 
The fixed points of Eqs. (40) that contain the additional 0{^/^) terms are q^ = {—p^,(f)s) 



and p^ = (pg, 0c), where p^ = ^/^3h^/T^^l^ /2. For small positive ^, a linear analysis of these 
fixed points reveals that q^ is still a saddle but that p^ is a sink, as shown in Fig. |4Fb). The 



fixed points of the full equations (39) near I = V are the same saddle and sink, located at 
(/ = J~, = (ps) and (/ = J"*", = 0c), respectively, where I^ = F ±: V^Pi- 



The scaled equations (41) provide an estimate for the basin of attraction of the sink, 



which is the area confined within the homoclinic orbit connecting the saddle go to itself, 
shown in Fig. |4|^a). Recall that the dynamics on the unperturbed annulus ^ is composed 
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(a) e = 



(b) C = 1 



FIG. 4. (a) Numerical phase portraits of Eqs. (40) with ^ = [or equivalently, Eqs. (41)], showing 



a saddle and a center, (b) Numerical phase portraits of Eqs. (40) with .^ = 1, showing that the 
saddle remains a saddle but that the center becomes a sink, with their p coordinates shifted slightly 
down and up, respectively. The parameters are 5 = 2,VL = 400, h = l,b = 0.2649, 7 = hb5. 



of simple one-dimensional flows, which on resonance turn into a circle of fixed points. Upon 
adding the small perturbation, two of these fixed points persist in an interval of length 
vr, and the phase space contains two-dimensional flows. Of particular interest is the basin 
of attraction of the sink, because a homoclinic orbit to a flxed point of this type offers a 
mechanism for producing chaotic motion. This mechanism, which results from the existence 
of a homoclinic trajectory to a saddle-focus flxed point, was described by Silnikov |35] . 
Obtaining an estimate for the basin of attraction of the sink, allows us to pick-out the 
trajectories satisfying Silnikov's theorem, which we do in the following section. 

VII. A HOMOCLINIC CONNECTION TO THE SINK p^ 

We are flnally in a position to show the existence of an orbit homoclinic to the sink 
p^. Note that for a particular set of parameters the existence of such an orbit implies the 
existence of another symmetric orbit due to the invariance {x,y) — )■ (— x, —y). To achieve 
this, we flrst show that there exists a homoclinic orbit that approaches p^ asymptotically 
backward in time, and approaches the perturbed annulus ^^ asymptotically forward in 
time. We then estimate the conditions under which the perturbed counterpart of the point, 
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which is reached on ^ forward in time in the unperturbed system, hes within the basin of 
attraction of the sink p^ on ^^. This gives us an estimate for the possibihty of obtaining a 
Silnikov orbit that connects the sink back to itself. 

The first step is done by finding the conditions for which the Melnikov function M{r', 0o = 
(pc + A0/2) has simple zeros. We substitute 0o = 0c + A0/2 into the first term in Eq. (36), 
and recall that sin 20c = —b, to get 



sin 200 sin A0 = - Vl -6^(1 - cos 2A0) - b sin 2A0 



(43) 



By substituting (43) into the Melnikov function (36) and equating it to zero we obtain the 
equation 



sr 



Vl-62(l - COS2A0) - 6sin2A0l + 26(3rAxi + 3Afi - 2Aa) = 



(44) 



from which we extract an explicit expression for the condition on b, ensuring the existence 
of an orbit that asymptotes to p^ backwards in time, and to ^^ forward in time, 

1 - cos 2A0 



1^1 



(45) 



(gl^Aa + sin 2A0 - 2Axi - ^) + (1 - cos 2A0)2 

Next, we wish to find an approximate condition, ensuring that this orbit approaches p^ 
as T — 7- oo. To do so we find the condition for which the unperturbed heteroclinic orbit, 
which asymptotes to po as T — )■ — oo, returns back to a point on the circle of fixed points 
that is inside the homoclinic separatrix loop connecting the saddle go to itself [3lj. Such an 
orbit is shown in Fig. [5j This condition is formulated in terms of the difference A0 between 
the asymptotic values of the angular variable as 



< 0c + A0 < 



(46) 



where 0m is the maximal value of on the homoclinic orbit, connecting the saddle go to 
itself. Since the Hamiltonian is conserved along an orbit, 0m satisfies the equation 



= e^(0, 

= SI'M 



m)-^(O,0.) 

-Vl- 62 + -COS 20, 



6 0m + 



TT 



arcsin b 



(47) 



whose roots are found numerically to obtain 



Eqs. (45) and (46) define conditions for the existence of orbits homoclinic to the sink p^. 
We wish to relate these results to the actual physical parameters of the coupled resonators. 
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FIG. 5. The heteroclinic orbit given by Eqs. ( 15 ) with / = P', superimposed with the phase portrait 
of the unperturbed scaled system on ^f near resonance, given by Eqs. (41). The parameters are 



the same as in Fig. Ha), with ili = 21.32. For these parameters b = 0.2949 according to Eq. (45), 
so we fix /i = 1 and 7 = 5bh in Eqs. (41 ). This value of 6 sets (j){—oo) = (j)c = —0.1341, and as can 
be seen from the figure (ps < 0(oo) = (j)c + ^4> < 4>m,- 

Recall that 5 sets the value of the electrostatic coupling coefficient D = 2{5'^ — l)/(35^ ~ !)• 
The scaled frequency Vli is then given by Vli = (wi — a;2)/2e = (a/1 — D/2— a/1 — 3D/2)/2e, 
so by fixing e it is also determined by 5. The ratio h = '~f/h6, between the damping coefficient 
and the drive amplitude, has to be positive in order for the damping coefficient 7 to be 
positive and have the standard physical meaning of energy dissipation. The ratio b is positive 

if the inequality 

4 . . . 2A/i ,, , 

—Aa + sin 2A0 - 2Axi - ^ > (48) 

is satisfied. We plot the left-hand side of this inequality as a function of fi and 5 in Fig. [6ta), 
and find that it is positive ii 5 > 1.5. Consequently we plot the ratio b in Fig. [6Fb) for 



1.5 < 5 < 3. This value of b then determines the values of the fixed points of Eq. (41) 



which are shown in Fig. l7](a), along with (pc + Acp and (pm for a particular value of il. The 
parameter values for which these values satisfy the condition (46 ) are displayed in Fig. [71(b), 
which outlines the values of the electrostatic coupling and parametric driving frequency, for 
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FIG. 6. (Color online) (a) Contour plot of the left-hand side of the inequality (48). In the displayed 



range of il., for 5 > 1.5, this function is positive and the coefficient 7 represents energy dissipation. 
For fixed e = 0.01 and 1 < 6 < 3, the scaled frequency f^i reaches values of < f^i < 28. (b) The 



ratio b = ^/h6, given by Eq. (45), as a function of il and 6 {e = 0.01). Here 1.5 < 6 < 3 and b is 
positive. 



which orbits homoclinic to the sink p^ exist. We note that Silnikov orbits were also found 
in other two-mode parametrically driven systems [391 SOI US], however, slightly different 
equations were studied, resulting in different phase space dynamics for the unperturbed 
system as well as different perturbations. 



Finally, we wish to verify our calculations by a numerical solution of the ODEs ( 19 ). The 



difficulty in producing a Silnikov orbit in these equations is that the linearized growth rates 
of the saddle- focus fixed point — a saddle on the {B, 6) plane and a focus on the perturbed 
annulus ^^ — are 0(0 iii directions tangent to ^^, so the orbit has to spend a lot of 
time near ^^ in order to spiral around the saddle-focus. However, the linearized growth 
rates of this fixed point in directions transverse to ^g, are 0(1), so a small and inevitable 
numerical error would defied the orbit away from ^^. To avoid this problem we solve the 



ODEs (19) using a cutoff criterion. We initiate the numerical solution with 5^1, and the 
exact coordinates of the sink on ^g, {I = I~^ , (p = (pc)- The orbit initially flows away from 
^^ and later turns around and approaches it. If on its way back towards ^^, the orbit 



approaches it close enough to satisfy B < ,^/1000, we set dB/dT = dO/dT = in Eq. (19), 
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(a) 



(b) 



FIG. 7. (a) (Color online) The values of i 



's, Vc V^c 



+ A(/), and (pm as functions of (5, for VL = 1135.64. 



For 6 > 2.12 the condition (46) is satisfied and orbits homoclinic to the sink p^ exists, except when 



A(f) = 0. (b) Parameter values for which the condition (46) is satisfied are indicated in gray. The 



white line inside the gray area corresponds to A(/» = 0, where the theory does not apply. In both 
figures e = 0.01. 

thus restricting the motion to be tangent to ^^. This numerical scheme allows us to verify 
our predictions, because as shown in Fig. |8} only when the damping coefficient is equal to 
7 = hb6 (± ~ 0.1%), with b given by Eq. (45), is our cutoff criterion for eliminating the 



motion transverse to ^^ satisfied. Furthermore, as shown in Figs. |8](a) and (d), among 



the orbits that satisfy our cutoff criterion, only the ones that satisfy the condition (46) 
asymptote to the saddle-focus. 

Owing to a theorem of Silnikov [M^, the existence of orbits homoclinic to a saddle- focus 



fixed point in Eq. (19) implies that these equations contain chaotic motion in the sense of 



horseshoes in their dynamics. 



VIII. SUMMARY 



We have studied the origin of chaotic dynamics, and provided conditions for its existence, 
in a case of two parametrically-driven nonlinear resonators. This was achieved by applying a 
method of Kovacic and Wiggins on transformed amplitude equations that were derived from 
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(a) 



(b) 





(c) 



(d) 



FIG. 8. Results of our numerical scheme for ^ = 0.001 and the rest of the parameters as in Fig. [5| 
(a) An illustration of the Silnikov orbit that is obtained for 7 = hb5 (The ^ — t- limit of this orbit 
is shown in Fig. [5J). For (b) 7 = /i(6 - 0.0003)5, and (c) 7 = /i(6 + 0.0002)5, and we see that the 
orbit does not get close enough to ^^ in order to meet our cutoff criterion for being homoclinic to 
it. (d) For i7 = 800, we obtain an orbit that approaches ^^ by setting the appropriate value of 7 
{b = 0.6178), however, this orbit does not asymptote to the saddle-focus, in agreement with Fig. M 
(b). In this simulation, the orbit leaves ^^ through its boundary at / = Qi/{6 — 1) and eventually 
/ — )> and the motion dies out. 

the equations of motion, which model an actual experimental realization of coupled nanome- 
chanical resonators. We considered the amplitude of the drive and the damping to be small 
perturbations and obtained explicit expressions for orbits homoclinic to a two-dimensional 
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invariant annulus in the unperturbed equations. At resonance, we were able to calculate the 
Melnikov integral analytically, and provide a primary condition for having homoclinic orbits 
in the full, perturbed equations. By further studying the effects of perturbations on the 
invariant annulus near resonance, we found a secondary condition for the existence of orbits 
homoclinic to a fixed point of a saddle-focus type. We used a numerical scheme to verify 
our theoretical predictions. Such Silnikov homoclinic orbits give rise to a particular type of 
horseshoe chaos, which can be expected in the dynamics of the full system for parameter 
values in the vicinity of those presented here. 
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